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^ Abstract 

<N 

We obtain a five-step approximation to the quasiperiodic dynamic 
^_H scaling function for experimental Rayleigh-Benard convection data. 

^-H When errors are taken into account in the experiment, the f(a) spec- 

^^ trum of scalings is equivalent to just two of these five scales. To over- 

come this limitation, we develop a robust technique for extracting the 
scaling function from experimental data by reconstructing the dynam- 



o 



^^ ics of the experiment. 

O 1 Introduction 

c^ 

"^ The most complete invariant description of a chaotic dynamical system is 

the dynamical scaling function [1] which provides scale factors that allow 
the reconstruction of the dynamics and the direct evaluation of its ergodic 
measures. Given the scaling function one can calculate all other invariant 
quantities (long time averages) describing the dynamical system, such as the 
spectrum of singularities f(a) or the correlation functions. Thus, the scal- 
ing function is the quantity of choice for characterizing a particular chaotic 
dynamical system. Theoretically the scaling function has been computed 
for the period-doubling and quasiperiodic transitions to chaos. Because of 
the completeness of the scaling function, a comparison between theoretical 
and experimental scaling function is a much more rigorous test of universal- 
ity than is commonly demonstrated with the f(a) spectrum of singularities. 



Despite the importance of constructing the scaling function from experimen- 
tal data, few attempts have been made and these have not been convincing 
[2]. The difficulty lies in the sensitivity of the scaling function to varia- 
tions in parameters such as can come from experimental noise or drift and 
in the proper definition of the scaling function, none of which have been 
addressed in experimental data analysis. We explain here how to overcome 
these difficulties and give the first definitive comparison of experimental and 
theoretical scaling functions for the quasiperiodic transition to chaos. 

Why have previous attempts to extract the scaling function failed? There 
are basically two reasons: The first is that the universal scaling function 
has to be computed not from any map, but from the universal function that 
satisfies the Cvitanovic-Feigenbaum functional equation [3]. Experimentally 
the universal function is not available, but it can be approximated from its 
definition by considering not the map obtained from the data, but one of 
its iterates. The second reason is that the scaling function is very sensitive 
to variations in the control parameters and these variations can never be 
completely eliminated. Numerical studies with the sine circle map show that 
large variations can be expected in the scaling function if the parameters 
are not exactly tuned to the quasiperiodic state. Here we will show how 
to approximate the behavior of the universal function by considering only 
orbit points close to the critical point, and how variation of parameters can 
be controlled by reconstructing the dynamics of the system. 

The circle map scaling function is a generalization of Shenker's contrac- 
tion rate a [4] to all points in the neighborhood of the infiection point of the 
circle map. Shenker's a measures the rate of exponential contraction of the 
close return distances of the infiection point of the circle map. We focus on 
the scaling function for circle maps, as it is a common case in physical sys- 
tems, arising generically when two oscillators are non-linearly coupled. We 
use the sine circle map for numerical comparisons and obtain experimental 
data from a hydrodynamical experiment: Rayleigh-Benard convection in a 
■^He-superfiuid-'^He mixture. This system, which closely approximates clas- 
sical thermal convection, has been extensively studied in the quasiperiodic 
regime [5, 6, 7]. 

The universality theory for circle maps is of wide interest because it 
occurs whenever two oscillators are nonlinearly coupled (the frequency of 
oscillation depends on the amplitude). If the coupling is strong, the sys- 
tem will go chaotic, but for any coupling there will be mode-locking. This 
was first reported by Christian Huyghens in 1665 when he described how 
clocks set on a shelf would synchronize the motion of their pendula [8]. The 



y y 

y y 
y y 
y y 
y y 
y y 



y y y 
y y y 
y y y 
y y y 
y y y 
y y y 
y y y 
y y y 
y y y 
y y y 
y y y 
y y y 
y y y 



y y y 
y y y 
y y y 
y y y 
y y y 
y y y 
y y y 
y y y 
y y y 
y y y 
y y y 



y y y y y 
y y y y y 
y y y y y 
y y y y y 
y y y y y 
y y y y y 
y y y y y 
y y y y y 
y y y y y 
y y y y y 
y y y y y 
y y y y y 
y y y y y 



+ 8 



i\ \ \ 
\ lu 

n\ r 
/ / /^ 
x/ / _ 

/ / / /^I:/ / f / , , 1 
\ \ 1 / >/ \ \ \_ I I I 

1 1 / ^r"-// 1 / ^, / 



■: \ r u j:^ 
. . , , ^ 



^/ / 



/! / ^1/! / / 



\ 




Figure 1: Mode locking occurs for typical vector fields on a torus. If a harmonic 
oscillator with irrational frequency is slightly perturbed, the winding trajectory will 
form periodic orbits. 



phase space of the two oscillators is composed of their amplitudes and phases 
and is thus four dimensional. In general, the oscillators are dissipative and 
therefore the study of their long term behavior may be reduced to the study 
of limit cycle sets. In the simplest case one of the oscillators has a fixed 
frequency and is driving the other oscillator, and the phase space may be 
reduced to the motion on the two-dimensional surface of a torus. On this 
torus there is a vector field that determines the motion. If it where a har- 
monic oscillator, then the vector field would be a set of vectors all pointing in 
the same direction, as the vector field (a) in figure 1. The typical frequency 
of oscillation is an irrational number (the rationals have measure zero), and 
the system winds around the torus ergodically. As the orbit winds around, 
it comes arbitrarily close to any point it has already visited. If we add a 
perturbation to this oscillator (vector field (b) in the figure), as we do when 
we couple it nonlinearly to another oscillator, then whenever there is a close 
return of the orbit to previously visited points there is a chance that the 
orbit will close on itself and be a periodic orbit — the system would be 
mode-locked. Peixoto [9] proved that this is the general case by showing 
that a typical vector field on a torus is mode-locked. 

Here we show a non-trivial physical realization of the mode-locking phe- 
nomenon. Peixoto's theorem applies to a two-dimensional dynamical sys- 
tem, but our analysis of the experimental data will show that the phe- 
nomenon is also observed in a system that is described by a partial differen- 
tial equation (Navier-Stokes equation) and therefore potentially an infinite 
dimensional dynamical system. The experiment is not used as an analog 
computer to simulate equations that are known to have mode-locking. A 
detailed analysis of the parameter space of the experimental system shows 



that it is not globally equivalent to the sine circle map [5]. There are non- 
chaotic regions in the experiment that do not occur in a simple circle map. 
This is not an unusual situation for higher dimensional maps that have in- 
variant circles (see, for example Wang et al. [10]). 

This paper is organized by first presenting a review (section 2) of the 
major features of the dynamical scaling function, emphasizing connections 
to physically obtainable data sets (either by experiment or numerical simu- 
lation). Next, in section 3, a brief description of the experiment is given. In 
section 4, we present the results of our data analysis, illustrating potential 
difficulties using numerical simulation, but concentrating on obtaining a re- 
liable scaling function for experimental data. Our summary and conclusions 
are contained in section 5. 

2 Scaling functions 

Fractals are complicated sets to describe. As a consequence several possible 
descriptions have been proposed, with varying degrees of completeness. The 
coarsest description of a fractal is its fractal dimension. It gives an idea of 
how many small boxes of fixed size are needed to cover the set. The closer 
the fractal dimension is to the embedding dimension, the closer it appears to 
be a figure of non-zero measure. A more detailed description of the fractal 
is given by the spectrum of singularities, the f{a) curve, or equivalently, the 
generalized dimensions Dq. The f{a) curve generalizes the fractal dimen- 
sion by decomposing the fractal set into self-affine fractals (which are not 
multifractals) indexed by a and for each a gives its fractal dimension. Most 
fractals encountered in physics have this multitude of scales and a parabola- 
shaped f{a) curve. The f{a) curve has more information than the fractal 
dimension, as it describes the decomposition of the fractal set. Despite the 
infinite number of fractal dimensions it contains, it is still not possible to 
reconstruct the fractal from the f{a) curve — /(«) is not a complete de- 
scription of the fractal. One might argue that a complete description of the 
fractal is not desirable, because it would be too complicated and because 
in principle the dynamical system that generated the fractal already pro- 
vides a complete description. Such a complicated description would not be 
useful for comparing experiments to theories. Such an argument would be 
correct if there where no structure to a fractal, but there is. Fractals that 
occur in physical system are seldom arbitrary, and are usually described by 
a smooth presentation function, or equivalently a scaling function. It was 



Feigenbaum who observed that there is a hierarchical structure to the de- 
scriptions of a fractal that can be explored to create a function — the scaling 
function — which can be easily approximated by a simple function. By a 
simple function we mean a function that has a good approximation in terms 
of a basis of computable functions. For example, most f(a) curves can be 
very well approximated (to less than 1%) by a parabola, and therefore are 
well approximated by three numbers and the basis functions 1, x, and x'^. 
The scaling function would be of little practical value if it were not well 
approximated in a simple basis, step functions in this case. 

There are many routes that lead to an explanation of what a scaling 
function is and how to compute it. The shortest is by breaking away from 
the historical development and considering first the presentation function of 
a fractal. The presentation function is a simple chaotic dynamical system 
(hyperbolic, unlike the circle map) that generates the fractal and is closely 
related to the definition of fractals of Hutchinson [11] and the iterated dy- 
namical systems introduced by Barnsley and collaborators [12]. From the 
presentation function one can derive the scaling function, but we will not 
do it in the most elegant fashion, rather we will develop the formalism in a 
form that is directly applicable to the experimental data. 

In the upper part of figure 2 we have the successive steps of the con- 
struction similar to the middle third Cantor set. The construction is done 
in levels, each level being formed by a collection of segments. From one level 
to the next, each "parent" segment produces smaller "children" segments 
by removing the middle section. As the construction proceeds, the segments 
better approximate the Cantor set. In the figure not all the segments are 
the same size, some are larger and some are smaller, as is the case with 
multifractals. In the middle third Cantor set, the ratio between a segment 
and the one it was generated from is exactly 1/3, but in the case shown in 
the figure the ratios differ from 1/3. If we went through the last level of 
the construction and made a plot of the segment number and its ratio to its 
parent segment we would have a scaling function, as indicated in the figure. 
A function giving the ratios in the construction of a fractal is the basic idea 
for a scaling function. Much of the formalism that we will introduce is to be 
able to give precise names to every segments and to arrange the "lineage" 
of segments so that the children segments have the correct parent. If we do 
not take these precautions, the scaling function would be a "wild function", 
varying rapidly and not approximated easily by simple functions. 

To describe the formalism we will use a variation on the quadratic map 
that appears in the theory of period doubling. This is because the combi- 
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Figure 2: Construction of the steps of the scaling function from a Cantor set. From 
one level to the next in the construction of the Cantor set the covers are shrunk, 
each parent segment into two children segments. The shrinkage of the last level of 
the construction is plotted and by removing the gaps one has an approximation to 
the scaling function of the Cantor set. 



natorial manipulations are much simpler for this map than they are for the 
circle map. The scaling function will be described for a one dimensional 
map F as shown in figure 3. Drawn is the map 
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restricted to the unit interval. We will see that this map is also a presentation 
function. 

It has two branches separated by a gap: one over the left portion of the 
unit interval and one over the right. If we choose a point x at random in the 
unit interval and iterate it under the action of the map F, equation (1), it 
will hop between the branches and eventually get mapped to minus infinity. 
An orbit point is guaranteed to go to minus infinity if it lands in the gap. 
The hopping of the point defines the orbit of the initial point x: x \-^ xi \-^ 
X2 1-^ • • •. For each orbit of the map F we can associate a symbolic code. 
The code for this map is formed from Os and Is and is found from the orbit 
by associating a if Sf < 1/2 and a 1 if Sf > 1/2, with i = 0, 1, 2, . . .. 

Most initial points will end up in the gap region between the two branches. 
We then say that the orbit point has escaped the unit interval. The points 
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Figure 3: A Cantor set presentation function. The Cantor set is the set of all 
points that under iteration do not leave the interval [0, 1], This set can be found 
by backwards iterating the gap between the two branches of the map. The dotted 
lines can be used to find these backward images. At each step of the construction 
one is left with a set of segments that form a cover of the Cantor set. 



that do not escape form a Cantor set C (or Cantor dust) and remain trapped 
in the unit interval for all iterations. In the process of describing all the 
points that do not escape, the map F can be used as a presentation of the 
Cantor set C, and has been called a presentation function by Feigenbaum 
[13]. 

How does the map F "present" the Cantor set? The presentation is done 
in steps. First we determine the points that do not escape the unit interval 
in one iteration of the map. These are the points that are not part of the 
gap. These points determine two segments, which are an approximation to 
the Cantor set. In the next step we determine the points that do not escape 
in two iterations. These are the points that get mapped into the gap in one 
iteration, as in the next iteration they will escape; these points form the two 
segments Aq ' and A-^ ' at level 1 in figure 3. The processes can be continued 
for any number of iterations. If we observe carefully what is being done, we 
discover that at each step the pre-images of the gap (backward iterates) are 
being removed from the unit interval. As the map has two branches, every 
point in the gap has two pre-images, and therefore the whole gap has two 
pre-images in the form of two smaller gaps. To generate all the gaps in the 
Cantor set one just has to iterate the gap backwards. Each iteration of the 



gap defines a set of segments, witli tlie ra-tli iterate defining tlie segments 

(n) 

A^ ' at ievei n. For tliis map tliere wiii be 2" segments at ievei ra, witli tlie 
first few drawn in figure 3. As ra ^ oo tlie segments that remain for at least 
n iterates converge to the Cantor set C. 

The segments at one level form a cover for the Cantor set and it is 
from a cover that all the invariant information about the set is extracted 
(the cover generated from the backward iterates of the gap form a Markov 
partition for the map as a dynamical system). The segments {A^"'} at 
level n are a refinement of the cover formed by segments at level n — 1. 
From successive covers we can compute the trajectory scaling function, the 
spectrum of scalings f(a), and the generalized dimensions. 

To define the scaling function we must give labels (names) to the seg- 
ments. The labels are chosen so that the definition of the scaling function 
allows for simple approximations. As each segment is generated from an 
inverse image of the unit interval, we will consider the inverse of the pre- 
sentation function F. Because F does not have a unique inverse, we have 
to consider restrictions of F. Its restriction to the first half of the segment, 
from to 1/2, has a unique inverse, which we will call F^ , and its restric- 
tion to the second half, from 1/2 to 1, also has a unique inverse, which we 
will call F^ . For example, the segment labeled A'^'(0, 1) in figure 3 is 
formed from the inverse image of the unit interval by mapping A^^^, the 
unit interval, with F^ and then F^ , so that the segment 

A(2)(0,l) = F-i(Ffi(A(°) 

The mapping of the unit interval into a smaller interval is what determines 
its label. The sequence of the labels of the inverse maps is the label of the 
segment: 



A(-\e„e2,...,e^) = F-' oF-' o...F-Ua(' 



3) 

The scaling function is formed from a set of ratios of segments length. 
We use I • I around a segment A(")(e) to denote its size (length), and define 






|AW(61,62,...,' 



|A(-1)(62,...,' 



We can then arrange the ratios (T'"'(ei, 62, • • • , ^n) next to each other as piece- 
wise constant segments in increasing order of their binary label ei, 62, • • • , e™ 



so that the collection of steps scan the unit interval. As ra ^ oo this col- 
lection of steps will converge to the scaling function. In section 4 we will 
describe the limiting process in more detail, and give a precise definition on 
how to arrange the ratios. 

The construction we gave for the scaling function cannot be used for 
the circle map or the quadratic map (the map ax(l — x), a < 4) because 
neither is hyperbolic. The essential point of the construction of Feigenbaum 
and collaborators was to realize that there is a way of re- writing these maps 
so that they are effectively hyperbolic. In both cases a universal function 
is constructed, and from it the scaling function or a presentation function 
can be computed. The universal function can be computed from the circle 
map /. Assuming that the map / has an infiection point at 0, Shenker [14] 
observed that if we compose / with itself a Qn Fibonacci number of times 
and choose a suitable value for a we can approach the universal function g 

a'^f^"{a-'^x)^ g{x) (3) 

as a; ^ and ra ^ oo. From this relation we discover that the universal 
function satisfies a functional equation 

Tg(x) = ag(ag(a~ xj) (4) 

(the usual functional equation uses the function ag). We interpret the a; ^ 
condition as stating that the universality results hold at the origin. The 
ra ^ 00 condition we interpret from the combinatorics of the circle map as 
meaning that we must iterate the map until the orbit lands close to the 
origin. For the computation of the scaling function we do not need, in 
principle, to restrict ourselves to the infiection point, as the scaling function 
is invariant under smooth diffeomorphism and also the action of the circle 
map itself. So if we use an image of the infiection point under the action 
of the map, we should be able to compute the scaling function. This would 
be correct if we could also take the ra ^ oo limit, but in practice data 
sets are limited. It is then no longer true that the scaling function can be 
computed at any image of the infiection point. A scaling function computed 
at an image converges more slowly than the scaling function computed at 
the origin. In section 4 we will give a procedure that effectively computes 
the scaling function from the universal function by using only iterates of the 
circle map. 



3 Experiment 

The experiment that provided the data for the analysis presented here has 
been well studied in the quasiperiodic regime using an apparatus that is de- 
scribed in detail elsewhere [5, 6, 7]. Here we give the essential features of the 
experimental data and discuss how the data are prepared for the calculation 
of the dynamical scaling function. The system is thermal convection in a 
■^He-superfluid-'^He mixture which approximates a classical convecting fluid 
with low Prandtl number. In some region of parameter space quasiperiodic, 
mode-locked, and chaotic states are observed. These states are the result 
of two internal oscillatory modes and not the result of external forcing. To 
study a quasiperiodic/mode-locking system one must vary two parameters 
independently. For this system, the two parameters are the temperature 
difference across the fluid layer and the mean temperature of the fluid. The 
range of rotation numbers that can be accessed by varying these parame- 
ters is from about 1/8 to 1/6. The canonical golden-mean rotation number, 
Pg = (v5 — l)/2 ~ 13/21 does not fall within this range but there are many 
rotation numbers with the proper golden mean tail (asymptotic series of Is 
in rational approximant series) that are in that range. For the purposes 
of testing universality of the quasiperiodic transition to chaos any of these 
rotation numbers is equivalent. The experimental data we use is centered 
around the golden-mean-tail irrational (3 — v5)/2. The strict golden mean 
has the advantage of making the "best" use of the data, as the data require- 
ments in terms of stability and precision do increase for rotation numbers 
other than the strict golden mean. 

The data are time sequences of temperature oscillations of the convective 
flow field measured at a local point in space. To the extent that this is a low 
dimensional dynamical system, measurement at a single point completely 
characterizes the state of the system. The data are used to reconstruct 
the phase-space dynamics using standard delay-coordinate embedding tech- 
niques [15, 16]. Poincare sections are produced by interpolating the inter- 
sections of the phase-space trajectories with a plane. For a quasiperiodic 
attractor, the section will fill up a smooth curve diffemorphic with a circle. 
In figure 4 we show a Poincare section for a state very near criticality (as de- 
fined in circle map descriptions of the quasiperiodic transition). Each point 
in the section has a time ordering according to its relative position in the 
time sequence and a space ordering that relates the nearest neighbor points 
along the Poincare section. The correspondence between the time and space 
ordering is determined by the rotation number and also sets the sequence 
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Figure 4: Poincare cross section from the reconstructed Bow. 

of close return times for the dynamics. Another effective action of the dy- 
namics on this curve is to partition the curve into segments that connect 
nearest-neighbor points. It is these segments that form the set which obeys 
multifractal scaling and from which we will construct the scaling function. 
Rather than considering the Poincare section, we can make a further sim- 
plification by constructing a one- dimensional mapping of arc length along 
the curve. Such a mapping, for the data in figure 4, is shown in figure 5, and 
is clearly a one- dimensional map very similar to a sine circle map. For more 
subcritical parameters, such maps constructed from data show all the simple 
features of a sine circle map [7]. Practical considerations that arise for ex- 
perimental data used to construct a multifractal description of the attractor 
(fractal dimension, f(a) spectrum, etc.) are the precision with which one 
can define a rotation number, the degree of random noise in the signal, and 
the stability of the state against drift in the parameters. In these experi- 
ments the signal to noise ratio was about 1000 : 1 and the rotation number 
could be determined to about 5 parts in a million. The most important 
factor which limited the data was drift in the operating parameter of the 
system. This caused changes in the rotation rate with time and, although 
very small, had extremely deleterious effects on the extraction of a scaling 
function for reasons discussed in the next section. In general the analysis 
for determining the scaling function is more demanding on the quality of 
the data than for averaged multifractal quantities because one is comparing 
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Figure 5: Return map reconstructed from figure 4 using arc lengtli. The dots are 
the data points and the solid curve the reconstructed dynamics as explained at the 
end of section 4. 

the ratios of individual segments as opposed to the averaging over many 
segments that defines the f(a) spectrum. In the next section we extract 
a scaling function for the quasiperiodic transition to chaos by reconciling 
the limitations of the experimental data with the correct definition of the 
scaling function. In particular, the proper limits must be observed. We 
demonstrate these methods on numerical sine circle map data to illustrate 
the pitfalls of scaling function analysis. 

4 The practice of circle map scaling functions 

In this section we will use the concepts developed earlier about scaling func- 
tions and adapt them to the requirements of circle maps and the realities of 
experimental data. We will have two types of difficulties. The simplest to 
overcome are those related to the combinatorics of the circle map, as it is 
not as simple as that of period-doubling maps. In the case of the circle map 
at a golden-mean-tail rotation number, not every segment gets sub-divided 
into two children segments; some do and others do not. The definition we 
adopt here for the scaling function matches that used by Feigenbaum in his 
presentation function article for circle maps [3], except that we use forward 
iterates. The other type of difficulties are associated with trying to deter- 



12 



mine the exact point where the data should be collected. We will conclude 
that it is not possible to collect the data reliably at the irrational rota- 
tion number and will therefore reconstruct the dynamics of the system as a 
function of one parameter in the vicinity of the irrational rotation number. 
In a numerical or laboratory experiment, owing to finite precision, the 
winding number is never an irrational, and the best that can be obtained 
is a rational approximant. In this case the map is locked into a periodic 
orbit and the range of parameters that have this frequency form a section 
of an Arnold tongue. The approximants are formed from truncations of the 
continued fraction expansion of the irrational winding number and form a 
series of fractions Pn/Qn- For example, if 

P= J = (ai,a2,a3,...) (5) 

ai -\ i 



(22 



as H 

then the approximants to the golden mean pg = (v5 — l)/2 = (1, 1,1,.. .) 
are: 

(i> = Y,(i,i> = ^,(i,i,i> = |--- (6) 

The numbers Qn are necessary to define the scaling function, and for the 
golden mean rotation number they are the Fibonacci numbers (Qo = 1, 
Qi = 2, Q2 = 3, Qn = Qn-i + Qn-2)- On the critical line, at the irrational 
winding number, one considers the first Qn points of the orbit. These points 

(n) 

delimit segments As along the circle, grouped in a series of levels, indexed 
by n. For the distorted loop from the experiment we calculate the separation 
of points using arc length along the curve. For the sine circle map the angular 
separation is used. The ratio of these segments is used to define the scaling 
function. In figure 6 the first 13 points of the irrational golden mean orbit 

in) 

are used to delimit the segments at level 3. The segments As are found by 

(n) 

iterating with the map the endpoints of segments Aq '. Because all points 
are obtained from the same orbit, the segments from previous levels can 
also be constructed from knowledge of the first 13 points, as indicated in 
figure 6. Notice that the ratio trick of reference [17] cannot be used, as all 
the different levels are computed at fixed parameters of the map, as required 
by the experiment. 

The scaling function is built from a series of piecewise constant steps of 
height (Ts placed in ascending order of the integer s (which indexes the 
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Figure 6: Arrangement of the segments used in constructing the scaling function. 

(n) 

The endpoints of the segments Aj. are the points indicated on the top line, and 
their ratios a^ are defined in equation 7. The levels n are indicated to the side, 

steps) and rescaled to span the unit interval. The scaling as' is given by 
the ratio between the size of the segment |As | and its parent segment: 

^(n)^J^ L (7) 

where Qn is the Fibonacci number of segments that are used at level n. The 
function 0(s, Q) is the parent index function which in the simple case of the 
golden-mean returns s — Q if s > Q and s otherwise. In figure 6, we have 
shown which segments to compare to compute the scalings as , as , and 

(2) . 

(Ts at different levels. Formulas for other rotation numbers are given in 
the appendix of reference [7]. The scaling function is a function of the unit 
interval into itself. In terms of the step function 9 (0(x) = 1 if a; > and 
otherwise): 

a(t) = J2ai:^)0(t-ty)0(t^J\!,-t) , (8) 

s 
(n) 

where the tl are where the discontinuities of the scaling function occur. 
The summation runs over a Qp Fibonacci number of them, and the value of 
p depends on the level of approximation to the universal scaling function. 
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For the simple case of a five-step scaling function they are: Iq' = 0, t-[^ = 

Qn-3/Qm't2 =Qn-2lQmt'i = Qra-1 /Qra, and ^4 = [Q n-3 + Q n-l) I Q n'l 

the general case can be worked out by expanding Mn a Fibonacci basis. 
To obtain a universal scaling function two limits must be considered. First 

in) 

n -^ 00, in which the scalings Gs go to their limiting values. In this 
limit more and more of the points of the irrational orbit are considered, and 
the segments As get closer and closer to the inflection point of the circle 
map. Second p —^ oo, in which the number of terms of the sum goes to 
infinity. The first limit (ra —^ oo) takes the scaling function to its universal 
form, whereas the second (p —^ oo) adds detail to the function. The limit 
towards detail cannot be taken before the limit towards universality. For 
experimental data the limit towards universality corresponds to considering 
a sequence of periodic orbits that approach the irrational one (tongue width 
going to zero), and the detailing limit corresponds to considering a larger 
number of levels. 

One practical consequence of the double limit is that one cannot use 
the smallest possible region determined by a periodic orbit as the segments 
for the scaling function. From numerical simulations we observe that the 
convergence is improved if we use only the first Qj points of a Qk > Qj 
periodic orbit in computing the ratios. For example, in an orbit of Qk = 
mil if we compute the scaling function with segments at level 8 the first 
step of the five-step scaling function is at 0.48, close to its limiting value 
of 0.468, whereas if we compute the scaling function with segments at level 
16 the same step is at 0.62, close to its trivial value. Notice that this 
implies that there are considerably fewer steps in the scaling function than 
would be expected from the experimental data set. In the scaling function 
we compute from experimental data we only consider the five-step scaling 
function. This is the largest number of scales we could extract given how 
closely we had approximated universality. As is illustrated in figure 7(a), 
the five-step theoretical scaling function is already a close approximation to 
the universal limiting function. 

As criticality is approached the scaling function goes from trivial behav- 
ior to behavior that characterizes golden mean criticality. In any experiment 
the rotation number is not exactly a golden-mean tail and the amount of 
data is not infinite, so for an experimentally determined scaling function 
the transition from the trivial case to the critical case is smooth instead of 
abrupt. This is analogous to the transition of the f(a) spectrum [18]. The 
effects of finite orbits and inaccuracies in the control parameters can be seen 
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Figure 7: The numerical and experimental scaling functions compared to the limit- 
ing universal scaling function. In (a) the five-step numerical approximation in gray 
is laid over the limiting scaling function in black. The two functions agree at the 
left end of the large steps. In (b) the five-step experimental approximation in gray 
is laid over the same limiting function. The error bars for the experimental curve 
can be found in table 1. 



in the scaling function. For a sub-critical orbit, if only the first few points 
of the orbit are used, then the scaling function resembles the critical sine 
circle map scaling function, but as more points of the orbit are considered 
the scaling function fiattens out, creating the apparent contradiction that 
as more data is considered, less "accuracy" is obtained. 

In figure 8 we have plotted several five-scale scaling functions for the 
sine circle map for a sub-critical value of the control parameter. The scaling 
functions differ by the number of points considered from the orbit. If only 
a few initial points are considered the scaling function resembles the critical 
one, but as more of the orbit is taken into account, the non criticality of the 
map becomes evident. Similar behavior is observed with the experimental 
data. The scaling function obtained resembles the theoretical curve as long 
as the orbits are short, but as longer orbits are considered the experimental 
scaling function differs from the predicted one. 

A similar phenomenon happens if we deviate from the superstable point 
along the critical line. In this case the scaling function is distorted away from 
the theoretical result, but still appears to be critical. In figure 9 we have 
plotted three scaling functions. One is for a rotation number smaller (below) 
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Figure 8: Several scaling functions computed at a small distance from the irrational 
winding number. The topmost curve is the theoretical curve and the next uses only 
the first few points of the orbit, whereas the bottommost uses around a quarter of 
the points. 
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Figure 9: Scaling functions along the away from the superstable point along the 
critical line. The plots show that large errors may occur when the system is not 
kept at the superstable point. 



the superstable point and the other is for a rotation number larger (above) 
than the superstable point. Notice that there can be large variations in the 
values of the scales. This points out the importance of being at exactly the 
superstable point in computing the scaling function. 

A straightforward application of the scaling function definition would 
consist of choosing a large-Q periodic orbit (a high order rational approx- 
imant to the golden mean) and using the smallest intervals at the largest 
level of an orbit to compute the scaling function. This disregards the order 
of the limits mentioned earlier. Such an approach was used previously, in 
combination with averaging, to reduce noise [2], and it fails because the con- 
vergence of ratios of intervals to their universal limit is nonuniform over the 
entire orbit and because averaging of noise does not improve uncertainty 
in the exact experimental winding number (in particular for drift in the 
control parameter, which we discuss later on). Explicit demonstration of 
the infiuence of noise on calculating the scaling function will be presented 
elsewhere [19]. 

Another problem with a straightforward application of the scaling func- 
tion definition is that segments approach their universal limit differently 
for different sections of the periodic orbit. Because the universality of the 
quasiperiodic scaling function stems from the self- similarity of any map un- 
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Figure 10: Segments around the infection point. This is an expanded view of the 
points around the segment from to 8. Tie extra points indicated in the top line 
are from the first 144 points of the orbit (the 89/144 approximant). Notice that 
the arrangement of the segments is the same as in figure 6. 



der composition and rescaling in the neighborhood of a cubic inflection point 
(the Cvitanovic-Feigenbaum functional equation), convergence is most rapid 
in the vicinity of this point. This translates in the computation of the scaling 
function into considering only segments around the inflection point, that is, 
not all the segments around the circle can be used as indicated in flgure 6. 
We accomplish this by considering the iterates of an initial segment not un- 
der the actual map _F, but of one of its iterates, F^^> = Fo- ■ -oF^ composed 
Qk times. In flgure 10 we have enlarged the region between the orbit points 
and 8 of flgure 6 and indicated the endpoints of the segments that would 
be used in computing the scaling function. The points of a 144/233 period 
orbit are plotted, but at most the flrst 89 points are used, which is three 
levels up from the bottom. This corrects the deflnition explained with the 
aid of flgure 6. 

In an experiment the circle map is not given analytically, and there is no 
direct way of determining the inflection point. We use the empirical criterion 
of Wang et al. [10]. It consists of observing that the inflection point, if it 
exists, is an endpoint of Aq ' which should be the largest segment at a 
given level. Although it is not essential that the segment containing the 
inflection point be used as Aq , the convergence is fastest if the segments 
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are organized around it. Plotting the return map or its derivative is not a 
method for determining the inflection point, as it may be absent in a map 
within an Arnold tongue [10, 20]. 

The flnal point that must be considered is experimental in nature and 
does not arise in numerical simulations. In order to obtain a reasonable 
scaling function one must be at a superstable orbit within a mode-locked 
tongue (in practice close to the middle of the tongue). In an experiment, the 
data set closest to the irrational winding number at the critical line is usually 
chosen. The closer the system is to the irrational winding number, the longer 
the orbit that can be obtained and the longer its control parameters must be 
kept locked within an Arnold tongue. But the width of a tongue decreases as 
the period of the orbit increases, and the longest orbit is obtained when the 
tongue width is below the resolution of the apparatus. In practice, for the 
narrowest tongues the system will jump between several mode-locked states 
as the control parameters are kept flxed within experimental resolution. 
To avoid this instability in the mode-locked state, orbits of shorter period 
should be considered. The experimentalist is then faced with the choice 
of either stable, short, and less converged orbits; or fluctuating, long, and 
better converged orbits. The reconstruction of the dynamics [21], described 
next, achieves an optimal compromise between these constraints and allows 
signiflcant improvements over the direct evaluation of the scaling ratios. 

To reconstruct the dynamics we determine two data sets that are close 
by in parameter space and on opposite sides of the golden-mean rotation 
number. We proceed to determine an interpolated map from each flnite set 
of experimental points using a least-squares cubic spline flt which is then 
iterated to determine its winding number. We then interpolate between the 
two maps to obtain the superstable point within one of the intermediate 
mode-locked tongue. (Because the two maps used for interpolation are close 
together, we use linear interpolation between their ordinates.) As a con- 
sistency check on our interpolation scheme we have computed the rate of 
contraction S of the Arnold tongues as the golden mean is better approxi- 
mated. It is measured to be 2.8, to be compared with the prediction of 2.83 
from the renormalization group for the circle map [4]. 

From the superstable interpolated map the flve-step scaling function is 
computed. With the interpolation method we can determine periods of 
lengths limited only by the computer, but we have been careful not to use 
periods that lead to average segment sizes that are smaller than the segments 
from the data. If this precaution is not taken the method will generate 
orbits whose universality class is dictated by the nature of the interpolating 
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Table 1: Five step scaling functions: limiting value, experimental value obtained 
with reconstruction, and directly without reconstruction. 



spline. The scaling function of the map reconstructed from the data is 
given in table 1 and plotted in figure 7(b). For comparison the theoretical 
scaling ratios for the universality class of the sine circle map, and the scaling 
function computed without the reconstruction process are also given. The 
theoretical scaling ratios were computed from a 832040/1346269 orbit of 
the sine circle map. The tabulated reconstructed scaling function is not the 
result of averaging over several data sets, but computed from a single long 
orbit. The errors are estimated based on several different rotation numbers 
with golden mean behavior. The effects of not being exactly at the irrational 
rotation number can be seen in the direct scaling function: the first three 
scales are in good agreement, but the final two, where small errors have 
accumulated, are not. 

5 Conclusions 

What have we gained in our analysis relative to, for example, computing 
an f(a) spectrum? First we know that the dynamics are correct because 
the construction of the scaling function requires constructing the symbolic 
dynamics of the map whereas the spectrum does not distinguish between 
fractal sets with the same statistics but different dynamics [22]. Second, we 
have extracted three scales beyond the f(a) spectrum [23]. Thus by extract- 
ing five scaling ratios that agree within 10% with the theoretical predictions 
we have made the most stringent test to date of quasiperiodic universality. 
In summary, it is possible to extract a scaling function from experimental 
data only if the orbit points around the infiection point are considered and 
if the parameters of the system are adjusted to be at the irrational winding 
number. Experimentally both constraints are interconnected and can be 
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resolved by reconstructing the dynamics. 
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